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ABSTRACT 


The present work is an attempt to carry out the nonlinear 
analysis of plates and linear analysis of shells. In both the 
cases , transverse shear deformation have been included. The 
nonl ineari ty is introduced by von-Korman strains . The response is 
obtained using first order shear deformation theory. Solutions 
have been obtained using finite element technique. A nine noded 
doubly curved isoparametric lagrangian finite element is used to 
d scretize the structure . The results are presented for linear 
bending , nonlinear bending, linear buckling and postbuckl ing of 
plates . However for shells only 1 inear bending and buckling 
analysis have been conducted. 

It is found that the first order theory ( FSDT ) may involve 
about E.5X difference with the higher order theory CHSDT) for a/h 
= 2 . After a/h = 10 there is hardly any significant difference 
between FSDT and HSDT which requires lot more computat ional 
effort . It is also found out that in case of buckling , FSDT theory 


gives results' which are close to 

Reddy *s 

HSDT 

results 

for 

all 

thickness ratios . This suggests 

that one 

need 

not go 

for 

the 


higher order theory in case of the buckling. 

To solve the nonlinear equations in the nonlinear analysi 
an incremental Newton— Raph son method has been developed which i 
more efficient than Pi card type of iterat ion. The theoretica 
results are compared with the experimental results and are foun 
to be in close agreement. The effect of different simply support 



boundary conditions has also been investigated. 

Separate programs written in FORTRAN are developed -for 
symmetric and unsymmetric laminates to reduce the computational 
effort . 



CHAPTER 1 


INTRODUCTION AND LITERATURE REVIEW 

1 . i Introductory Remarks s 

In recent years fibre reinforced composites , have found 
increasing use in many engineering structures and consumer 
appl ications, ranging from fuselage of an airplane to the frame of 
a tennis roquet. This is mainly due to two desirable features of 
reinf orced composites ie. , a high strength to weight ratio and 
flexibility in tailoring the anisotropic material properties by 
changing the fibre orientation, stacking sequence and thickness of 
laminae . The increasing use of composites as a structural element 
has generated considerable interest in the analysis of composite 
plates and shells. 

Composites can be considered as a material consisting of 
two or more chemically distinct constituents with a distinct 
interface between them. The properties of composites are 
noticeably different from constituent materials . Fibrous 
composites consist of essentially high strength fibres embedded in 
a matrix. A lamina is a thin sheet with fibers, oriented in 
certain directions, forming the main load bearing members. The 
matrix keeps fibres in position and transfers the stress to fibres 
through bonding. The lamina can be cons idered to behave as a two 
dimensional orthotropic material . 

The main advantage of composites lies in tailoring the 
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material properties to suit the applied loads. The disadvantage, 
from analytical point of view is the existence of coupling between 
bending and stretching. Due to this the governing partial 
differential equations for laminae become coupled in displacements 
and rotations . This boundary value problem is difficult to solve . 
Very few exact solutions can be obtained for limited boundary 
conditions . This is a marked departure from convent ional isotropic 
materials. 

Due to the high ratio of elastic modulus to shear modulus 
(of the order of 25 to 40, instead of 2.6 for a typi cal isot ropi c 
material ) , shear deformation can be significant in composite 
materials as compared to isotropic materials. The classical plate 
theory which does not take into account the transverse shear 
deformation, under predicts the deflections and over predicts the 
buckling loads and the natural frequencies , when the ratio of 
lateral dimension to the thickness is less than 50. 

For small deflections ie. , the deflections which are small 
compared to thickness, the linear theory gives sufficiently 
accurate results. Plates, subjected to lateral loads and/or in 
plane loads, fail at large deflections. For large deflections the 
lateral deflection is accompanied with mid plane stretching. This 
increases the load carrying capacity of the plate. Associated 
increase in load depends on the in plane boundary conditions 
imposed on the edges. This behavior is extremely important in 
structural weight minimization in aerospace structures. A plate 
starts to buckle when the critical load is reached. The maximum 
stress at the instant of buckling is much lesser than the yield 
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stress of the material. Hence the plate can take more load. As the 
in plane load is increased, the deflections increase at a faster 
rate. Now the deflections are no more comparable with the 
thickness, violating small deflection theory. Therefore large 
deflection theory should be considered to determine accurate load 
carrying capacity of laminated plates. 

The rapid development of computing technology has 
completely changed the out look of structural analyst. Of all 
approximate methods to analyze the structures, involving complex 
geometries and variety of boundary conditions, finite element 
method (FEM> is found to be most powerful . FEM, in structural 
analysis, is a technique that discretizes a structure in to a set 
of simple sub domains called finite elements . Over each element 
the solution of the governing differential equation is 
approximated by a 1 inear combination of undetermined parameters 
and pre selected algebraic polynomials. The undetermined 
parameters are determined such that the variational statement 
obtained from the governing differential equations is satisfied in 
each element. The satisfaction of the variational statement 
results in a set of algebraic relations among the undetermined 
parameters of an element. Since the domain is divided into a set 
of elements, they are put together into their original position, 
called assembly of elements, using continuity at element 
interface. The finite element method typically entails the 
solution of a very large number of equations for the nodal values 
of the function being sought. Thus the method requires so much 
computation that it is practical only if calculations are carried 
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on a computer. 

In the present work FEM is used to carry out the linear and 
non linear analysis of composite laminates including transverse 
shear deformation. 

1.2 Literature Review : 


The pertinent literature survey is presented here under the 
following three headings. 

Ci) Studies using small deflection theory of plates. 

(ii) Studies using large deflection theory of plates. 

< iii ) Studies devoted to small deflection theory of shells. 

1.2 CD Studies Using Small Deflection Theory of Plates s 


A number of theories have been proposed for shear 
deformable laminated plates. Stavsky M3 extended the classical 
plate theory for isotropic materials to laminated plates. This 
theory has been generalized to laminated anisotropic plates by 
Yang, Norris and Stavsky C23 and is called as YNS - theory. This 
is generalization of the Mindlin plate theory which includes 
transverse shear deformation and rotary inertia effects- In this 
theory the normal to the mid plane before deformation remains 
straight after deformation but not necessarily normal to the mid 
plane. To make shear strains zero on top and bottom faces of the 
plate, a shear correction factor is used which depends upon the 
ply lay-up and fibre orientation, etc. It has been shown that 



C3-53 YNS theory predicts the transverse deflections and natural 
frequencies for laminated composite plates with reasonable 
accuracy. Whitney and Pagano C63 presented some close form 
solutions for cross ply and angle ply rectangular plates under 
sinusoidal loading. Sinha and Rath £73 considered both vibration 
and buckling for the same type of plates. They investigated the 
influence of shear deformation , panel curvature and aspect ratio 
on bending stretching coupling. They used Wavier type solution to 
reduce the governing partial differential equations to algebraic 
equations. 

The finite element analysis of laminated composite plates 
began with Pryor and Barker C83 who employed twenty eight degree 
of freedom <dof ) shear deformable rectangular element. Pica C93 has 
compared the results of linear, serendipity, lag rang ian and 
heterosis elements with analytical thin plate solutions for 
square, ci rcular and elliptical isotropic plates based on YNS 
theory. He suggested the use of heterosis elements for rectangular 
and curved element meshes. Reddy C10U developed C° penalty plate 
bending element based on YIMS theory. The element consists of five 
dof at each node- three displacements and two rotations at the mid 
plane of the plate. Reddy C113 presented comparison of the finite 
element solutions using penalty method and exact solutions for 
anisotropic plates. He used reduced integration to correct the 
shear stiffness. He concluded that higher order elements are not 
affected by locking. He showed that nine noded lagrangian element 
is not much affected by the integration scheme but eight noded 
element is profoundly affected. The consideration of symmetries in 
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unsymmet ri cal laminated composite plates was investigated by Reddy 
C123. He presented the boundary conditions to be used in -full and 
quarter plate models of cross ply and angle ply laminated plates. 
Hass and Lee C13ZI used a nine noded degenerate shell element based 
on modified Hellinger Reissiner principle. Functions were assumed 
independently for in plane and transverse shear strains. 

The finite element formulation for the isotropi c plate 
buckling, using classical theory, is available in standard texts 
C14, 153 . Gaj bi r Singh and Sadasiva Rao C223 used an eight noded 
Mindl in element for the stabi 1 ity analysis of thick angle ply 
laminates. The effect of fibre or ientation , material properties, 
layering and boundary conditions on buckling loads were studied. 

Further development in laminate analysis have been higher 
order theories. Phan and Reddy C16D presented a higher order 
theory where in plane displacements are expanded as cubic function 
of the thickness coordinate while the transverse deflection is 
kept constant over the thickness of the plate. They presented 
closed form solutions for the bending of symmetric cross-ply 
laminates. Recently Phan and Reddy C17U using the higher order 
theory presented closed form solutions for the stability of 
isotropic, orthotropic and laminated plates. 

1.2 CiD Studies Using Large Deformation Theory of Plates i 


Much of research in the laminated composite plates is 
limited to linear problems. This is due to the complexity of non 
linear partial differential equations associated with large 
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deflection theory of composite plates. Approximate solutions for 
composite plates using large deflection theory based on Von Korman 
theory of isotropic plates have been attempted by C23-243. Chandra 
and Raju C233 and Chia and Prabhakara C243 employed the gelerkian 
method to reduce the governing di f f e rent ial equations to an 
ordinary differential equation, in time for the mode shape. The 
perturbation technique was used to solve the resulting equations. 
They analyzed bending and free vibration of composite laminates. 

Using FEM the total stiffness matrix can be written as the 
sum of linear stiffness matrix and non linear stiffness matrix. 
The non 1 inear stiffness matrix is called as incremental stiffness 
matrix. Two notations have emerged for the incremental matrices 
formulated in total lagrangian system. In total lag rang ian system 
al 1 the quantities of interest are expressed as functions of 
rectangular co— ordinates of the body in its initial or undeformed 
state. The first notation is due to Mallet and Marcal C251 and 
will be referred to as END notation. In this type of formulation, 
the expressions for total energy, and equilibrium are expressed as 

T T 

n = {«»} ( i L K ]„ + i [ N1] + 4 [ N2 ] ]{q} - {q} {«} 

( C K + ii l m 3 + is ["*]] { q } = { R } 

(-0 

Where K is the linear stiffness matrix and elements of Ml 
and N2 depend linearly and quadratically on the transverse 
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deflection u^. 

The incremental equilibrium equations can be written as 

[ [ « 3 t - [ N1 ] + [ «2 ] ] {•*<)} = {ar} 0-2) 

The second notation is due to Zienkiwicz and Nayak C263 

referred to as CBU notation In this type of formulation the 

equilibrium equation written as 

( [ «],. * [ K]_ + C ><L] {-} - {«} 0-3) 

where K is the linear stiffness K is unsymmetric stiffness 
matrix and K is the initial stress matrix, due to large 

Cf 

deformation theory. 

The incremental equilibrium equations can be written as 

[ K ]„ {q} = {arJ 04) 

where is the tangent stiffness matrix. The various matrices , 

K , K and K are derived for first order shear deformation 

ML Cf T 

theory in chapter 3. 

The correlation between the two notations was presented by 
Wood C273. 


Reddy and Chao C283 used mixed finite element for non 
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linear bending of composite plates in total Lagrangian system. 
Putcha and Reddy C293 used mixed finite element for the non 1 inear 
analysis of laminated plates based on higher order deformation 
theory. They used Picard type of iteration to solve the non linear 
algebraic equations with small load increments. 

1.2 Ciiil Studies Devoted bo Small Deflection Theory of Shells ; 


The analysis of laminated composite shells began with 
Ambartsumyan C303 . He assumed that individual layers are oriented 
so that principal axes of material symmetry are coincident with 
the principal coordinates of shell reference. Dong, Pister and 
Taylor C31U formulated a theory for laminated shells which is an 
extension of the developed by Stavsky C13. His work can be called 
as Donne 1 *s shallow shell theory. Cheng and Ho L3Z1 analyzed 
laminated anisotropic shells using Flugge's shell theory C331. The 
theory which considers shear deformation and thermal expansion 
through the thickness was considered by Zukas and Vinson C34U. 

The finite element analysis of composite shells was 
attempted by Pandya and Natarajan C353 and Shivakumar and 
Krishnamurthy C363 and Rao K.P. C373 based on Kirchoff's theory. 
Pandya and Natarajan C353 used a three noded on# dimensional 
curved finite element having four degrees freedom at each node for 
the analysis of axisymmetric fibre reinforced laminated shells, 
subjected to axisymmetric loading. The results for a cylindrical 
tube subjected to internal pressure and axial tension, are 
compared with the analytical solution. Shivakumar 


and 
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Krishnamurthy C363 used eight dof and sixteen dof ring element ■for 
free vibration analysis of laminated cyl indri cal and conical 
shells. They showed that for the same number of unknowns results 
with sixteen dof ring element are in good agreement with the 
experimental results than eight dof element. 

Reddy C381 developed shear deformable f ini te element for 
composite shel Is based on Sander * s theory C19H. He compared the 
approximate shell theories of Love , Sander, Donnel, Merely and 
Loo. He concluded that the effect of shear def ormat ion on the non 
dimensional deflections is more pronounced with the decrease in 
layers . He showed that there is a small difference in results 
obtained with different shell theories for shallow shells 
(R/h > 100 ). But for R/h = 10, the dif f erence in the results 
obtained with different shell theories is noticeable. Using Von 
K orman strains and Sander’s theory, and incorporating transverse 
shear deformation, Reddy et.al.C40H presented geometric non-linear 
analysis of doubly curved shells. 

1.3 Aim and Scope o t Present, Work i 


The literature survey presented reveals the necessity of 
nonlinearity and shear deformation in the composite plates. There 
are some higher order theories proposed till now to improve the 
solution accuracy for the linear problems. These are 
computationally expensive and involve large number of nodal degree 
of freedom. The analytical methods can be applied to linear 
problems with limited boundary conditions. But they are not 
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tractable for nonlinear problems. 

The scope of the present study can be summarised below 

In Ref . Him it is stated that higher order elements are 
not much effected by locking and nine noded lagrangian element is 
superior than eight noded element . Therefore the present study is 
an attempt to study the analysis of plates and shells using nine 
noded shear def ormable element . The displacement functions requi re 
C° continuity which leads to easy formulation for thick and thin 
plates. 

To compare the numerical reliability and numerical accuracy 
of the first order shear deformation theory CFSDT1 compared to 
higher order theory CHSDTI3 for linear bending and buckling 
problems. The FSDT is computationally inexpensive and easy to 
implement compared to HSDT . 

To check the validity of the Von-Korman nonlinear theory, 
the theoretical results are compared with the experimental 
results . 

A picard type of iteration is used in most of the works in 
solving the nonlinear equations C28,293. This type of iteration is 
computationally expensive. Therefore incremental Newton-Raphson 
method has been developed. 

In the nonlinear analysis the effect of different in plane 
boundary conditions is significant. The effect is more pronounced 
with bending- st ret ching coupl ing . Hence the nonlinear analysis of 
two layer cross-ply and two layer angle-ply with S1,S2 and S3 
boundary conditions given in chapter 2 is studied. 
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1.4 Layout of the Thesis s 


The thesis has been divided in to five chapters . The first 
chapter is an introduction to the problem with the pertinent 
literature. 

The second chapter deals with the theoretical formulation 
and der ivat ion of the expressions which can be readily used in the 
finite element formulation. THe finite element formulation and the 
solution procedure is presented in chapter 3. Chapter 4 deals with 
the interpretation of the results obtained and is followed by 
conclusions and scope of f urther work in chapter 5. 

The thesis has Appendix-I and Appendix— II. The figures and 
tables are given in the chapters that follow as and when the need 


arises. 



CHAPTER 2 


THEORETICAL FORMULATION 


2.1 Introduction : 

This chapter presents a detailed formulation of shear 
deformation theory of shells based on Sanders theory, including 
geometric nonlinearity. The geometric nonlinearity is introduced 
by von-Korman strains. The nonlinear partial differential 
equations for plates and shells cannot be solved analytically for 
complex geometries and various boundary conditions. FEM is the 
most effective numerical method to solve such problems. 

FEM is a piece wise application of variational method. In 
variational method the governing differential is recast in an 
equivalent integral form and a variational formulation is derived. 
This results weaker continuity requirements on the solution in the 
variational problem than in the original equation. 

In solid mechanics problems the variational form of the 
governing equations represents the total potential energy of the 
system. The solution of the differential equation is found by 
minimizing the variational form with respect to unknown variables. 
For monotonic convergence the primary variable should be 
continuous up to one order less than highest order, occurring in 
the variational functional. 

The classical plate requires cubic or higher order 
interpolation functions for approximating transverse deflection. 
But this first order shear deformation theory allows us to use 
linear or higher order inter polation functions for all 
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displacements . This type of formulation is easy to implement and 
applicable for both thin and thick laminates . However, shear 
deformable elements suffer from locking phenomenon. It is 
necessary that the shear strain fields, as the plate thickness 
reduces, should satisfy the Kirchoff *s condition ie. , 



= ^ 

+ w 

> 0 

xz 

X 

y X 


and y 


+ w 

> 0 

yz 

y 

tY 


Conventional methods 

using 

exact 

integration 


consistently, introducing spurious constraints that lead to 
locking and severe stress oscillations. Therefore reduced 
integration technique is used to circumvent the locking. 


2.2 Assumptions made in the Theory : 

Theoretical development in the present work is based upon 
the following assumptions. 

( i ) The laminate is constructed of an arbitrary number of 
layers of orthotropic sheets bonded together. However the 
orthotropic axis of material symmetry of an individual 
layer need not coincide with the laminated coordinate 
system. 

(ii) The strains are small. 

(iii) The transverse normal stress and strain are negligible. 

(iv) The line normal to the mid plane after deformation remains 
straight but not necessarily normal to the mid plane. 

(vJ There are no body forces. 

(vi) Stress-strain relations are linear . 



2. 3 Shell Geometry : 


-D 


Fig. 2. 1 shows a differential element of a shell. Here 
< £ _j * f £ ) denote curvilinear coordinates such that £ ^ and 
$ gCurves are lines of curvature on the mid surface C = 0. The £ 

curves are straight lines perpendicular to surface £ = O. The 

lines of principal radii of curvature coincide with the coordinate 
1 ines and the values are denoted by and R £ . X^ and X g denote 
cartesian reference frame. 

The 1 inear displacement field is given by, 

u C £ ^ £ z r £ ) = u lj ( « , ? 2 , 0 > + £ ^ 
v ( £ r £ rT £ ) = u 2 ( ? ^ £ zr 0 )+ £ (2.1) 

w ( £ „ ? 2 , C ) = u 3 ( e.j, ^ 2 , 0 ) 


Where u, v, w are the displacements along (? ^ , £ 2 , C ) coordinates, 
u , Ug, Ug are the displacements of the corresponding point on the 
middle surface and <p ^ and <j& 2 are the rotations at £ = O about Cg 
and £ ^axis respectively. 

The strain displacement relations including geometric 
nonlinearity can be written as 


‘1 

‘Z 

S 4 

e 5 

S 6 


o 
= 1 
a 
E 2 
o 
= 4 
o 
=5 
o 
“6 


£ 

£ 


+ £ K 


( 2 . 2 ) 



O Cl 

Where «. and k. are the mid plane strains and curvatures 
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respectively. They can be written as C403 
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& X j R 


( a _ 2 
3 1 




k. = 


S X. 


_o _ ^ U 2 
2 <9 x„ 


u. 


1 

2 


u 3 ' 2 
S X^ 


a 4>„ 


= 


a x. 




* u. 


* X, 


£ li. 




a u. 


$ X . 


}( 


<9 u„ 


0 X, 


s $ s 4> 2 

k 6 = — + — 


d X, 


# X 


( 


£ II, 


<3? x , 


a u 


a x 


2 ) 


& u_ 

o , 3 

*4 = + 

* x„ 


R, 


<9 u„ 

° _ j j. 3 

= = ^ + 

* X, 


u. 


.(2.3a) 


Where 


C o " 


1 

R, 


) 


(2. 3b) 


The term C is due to Sander's theory which accounts for the 
o 

condition of zero strain for rigid body motion. 
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2.4 Stress -Strain Relations for a Lamina * 

First order shear deformation theory used in the present 
work is incapable of satisfying the condition of zero transverse 
strain on top and bottom surfaces. Correction factors K^and 
have been used to account for this in stress strain relations . 

Based on the assumption that the transverse stress is 

negligible, the stress strain relations for an orthotropic lamina 
k oriented arbitrarily to the laminate axis are given by 



Eqn. (2-4) can be written in more compact notation as 

Hr [*].{•} 

Expressions for the elements of CQ3 matrix in terms of engineering 
constants for lamina for a given orientation a are given in 


Appendix— I. 
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2.0 Stress-Resultants at the Mid Plane of the Laminate s 


The stress resultants , N 2 ” N 6 T , M g „ and Q g of 

the laminate with «l number of layers is given by 


( U^, N gr N 6 ) 


HL 

E 

L-i 


K+i 




, tf 2 , & 6 > d K 


K+i 


1' 2 



Mia 

r 

n, > = 

6 

E 

Lsi 

| ( 1 ' & Z r *6 * ** 


Mia 

r 

) 

E 

( *5 } d K 


la* 4. 

£ 

N K 

resultants 

are 

as shown in Fig .2.2. 


(2-5) 


Substituting relations (2.2) in Eqn.2.4 the following 
expressions for stresses are obtained 




® jj <| 

°12 

°16 ‘ 


■-v 


' Q 11 

Q 12 

°16 ' 


f k 1 
1 

“ ^2 

k. 5S 


Q__ 

22 

°26 

- 

o 

*2 



Q 22 

*26 

K ■ 

k 2 ' 

- & b - 


Sym 


®66 - 


o 

^6 ^ 


Sym 


°66 - 


k° 

k 6) 



( 2 . 6 ) 
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-°5 


> = 


K f Q li 

1 44 


K 1 K 2 S 45 
4 a 55J 




Rewriting the stress resultants (2.5) in terms of mid plane 
strains and curvatures after integrating through each layer 


r N 1 ' 



r i 


r 1 

' 


1 

1 

N 2 



GA 


! GB 



o 



L J 


L J 



^2 

N 6 

p 

m 

• aw* 

a 

a a a 

m m 

m m 

H 

o 

► 





a 

a 

■ i 



4 

M 2 



Sym 

a 

a 

GD 



4 

‘ W 6- 




a 

a 




• k 6 ■ 


Cl ■ 


QS 


15 ) 


(2.7) 


Where GA, GB, GD and GS are extensional , bending - 
extensional coupling, bending, and shear stiffness matrices 
respectively. They are given by. 


f[GA], [GB], [0D]1-E J [Q] d( 

i = j = 1,3,6 (2.8) 


[ GS] 


£ 31 * 

E J [3] K 

1 ** 4 . £ 

K 

i = j = 4,5 


(2.9) 
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2.6 Boundary conditions s 


Since displacement based formulation is used in the present 
study only displacement boundary condition has to be specified. 
One can develop Navier solutions for different types of laminates 
with three types of simply supported boundary conditions . They are 
referred in this thesis as S1 r SS and S3 respect i vely . 

The functions which satisfy SI boundary conditions are 


u 1 (x 1 

T Xg ) 

II 

E 

m 

U cos 

mri 

otx . 
1 

SIM 

/5x 2 

U 2 <X 1 

r Xg 5 

= E 

r*=4 

2 

m*4 

V SZN 

mri 

OiKi 

COS 

/?x 2 

U 3 (X 1 

»x 2 > 

= e 

E 

4 

y sim 

mr» 

otx 1 

SIM 

ftx z 

*1 <X 1 

»« 2 > 

= E 

r»=4 

E 

W*= 4 

X cos 

mn 

SK, 

1 

SIM 


*2 <X 1 

* X g* 

= E 

nasi. 

2 

m=4 

Y SIH 

mri 

1 

COS 

ftx z 

N^o, 

x 2 ) = Nl 

( a , x 2 > 

= N 2 

(x^o) 

" N 2 ( *1 

,b) 

M^o, 

x 2 > = H, 

(a,x 2 > 

= M 2 

( X ^ y O 1 

- Vi 

yb) 


( 2 . 10 ) 


Where « 
shell . 


= — and ft = — t a and b are plane form dimensions of the 
a b 

The expressions which satisfy S2 boundary conditions are 


u 1 (x 1 

,Xg) 

= 2 

2 

U SIM 

mr» 

dX A 

1 

cos 

/5x 

u 2 <x 1 

,x 2 ) 

= 2 
r»*4 

2 

r«=4 

V cos 

a tx H 
1 

SIM 

/?x, 


N 6 (o,Xg) = N 6 (a»Xg) - N^x^yO) - N^tx^b) = 0 

M 1 (o,Xg) = (a,Xg) = M 2 <x 1 ,o) * Mgtx^b) = 0 
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The remaining expressions are same as given in Eqn.2.8 

The boundary conditions S1„ S2 and S3 for full shell 
quarter shell for m * n = 1 are as shown in Fig .2.3 

For clamped case al 1 the displacements and rotations 
restrained at the edge ie . , 

U 1 = u 2 = U 3 = ^ 1 = ^2 = 


and 


are 


0 


(2.13.) 








CHAPTER 3 


FINITE ELEMENT FORMULATION AND SOLUTION PROCEDURE 

3. 1 Introduction : 


In this chapter a displacement finite element formulation , 
based on total lagrangian for nonlinear analysis is presented. The 
continuum is divided in to number of elements and over each 
element a set of functions are chosen to define the state of 
displacements. The total potential energy is written in terms of 
displacements and using the principle of stationary of total 
potential for an element, we obtain equations which are nonlinear. 
These equations are solved by Newton-Raphson method. 

3. 3 El emani s 


The element used in the present work is a shear deformable, 
doubly curved shell element whose proj ect ion is an isoparametric 
rectangular element on plane. Over the typical shell 

element A , the displacements u^, u 2 , Ug, ^ and are 

interpolated by the expressions of the form 

NN 

-L - 2 N , ‘ v v < 


i 


1,3 


(3.1a) 
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MN 

e 





*<. ‘ 2 N „ ( 'l' *L 


L ss 

1,2 

(3.1b> 


j=1 





Where mm is the 

number of nodes (9) , N 

Jt 

is 

the 

shape 

functions 


e 

e 




associated with 

node j and u , 4> J 

i i 

are the 

unknown nodal 

variables. The 

shape functions are given 

in 

C143 




3. 3 Total L&grangian Virtual Work Equation s 


The basis of the present formulation is the virtual work 
equation for a continuum written in total lagrangian coordinate 
system under the assumption of small strains. 


I d «s. T <*. dv 


1 p du T q. dv 

* r 

du. T R. dA 

J L 1 


J t ) 

J 

J } 


* t= 1 ,5 

j-1,3 (3.2) 


Where v is the undeformed volume, p is the mass density, q is the 

v 

body force per unit mass and R. is the surface tractions per unit 
area, acting over the undeformed area. Here symbol d is used to 
denote first variation. 

The internal virtual work d 14 e over a element is given by, 

T 

d e. &. dv (3* 3a) 
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hs£ 


d W 


I I -V 


^ d < dA 


( 3 . 3b ) 


2 


Substituting the strain relations Eqn. (2.2) in to 
Eqn. (3.3b) r we obtain. 


K/2 


dW 


= I I [ “< -1 - C > 


<3' 1 + d( «g + ( kg ) ffg + 


-h^z 


d ( -c^ + £ k*J > <y_ + d o' + d <y_ d C dx - dx_ 

6 66 4 4 5 5 1 2 


(3.4) 


Where h is the thickness of the laminate. After integrating 
through the thickness of the laminate, Eqn. 3. 4 can be written as 


dW. 


I 


1 d + M 1 dk® + n 2 d ^2 + M 2 dk 2 + 


N 6 d m l + M 6 dk l * Q 1 d *4 + ®2 d *5 ] dX 1 dX J 


(3.5) 


Defining the mid plane strains, bending strains and shear strains 
in a matrix form as. 
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W 

W 


o o 
& 2 * & & 


{ *°- 
{ k °i’ k Z’ k 6 } 

{ *V “5 } 


(3.6) 


The mid plane strain « can be written as sum of linear 
strain and nonlinear strain. The linear strain is due to membrane 
strain and the nonlinear strain is due to bending deformation u„ 


only. 


Where 


{-r 




(3.7a> 


{-):■ 


'^1 + ^ ■ 


r 1 

. 2 

3) 



— 

— 

ix, R. 

1 1 


2 

»x 1 J 

<9u_ u„ 

r i° 

1 

au 3l Z 

^ mmm ^ 

r ? 1*1 = " 

— 

— 

d Kg Rg 

l Jml 

2 

J 

— + 


r ^3 i r ^ 3 ] 

. ax 2 * X 1 • 


l J l *x 2 J 


C3.7b) 


3.3.1 Expression For s 


€ 


- b; 


The membrane strain can be written in a matrix form in 

* L* 
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terms of unknown nodal displacements . Using the Eqn. (3. 1 ) in 
Eqn. (3.7b) one gets 




(3.8a) 


Here [ B is a matrix of size 3x45 and is given below 


[ B 3! - [ B, 


01 T B 02 t 


B 09 Ji <3 - 8b> 


B_ are the sub matrices and is given by, 
Oj 


[»]*. = 


«*N 

j 

Sx , 


&M 
J 


J 

£ 

dN 1 

J 

Sx d 


N e 


N. 

} 

R, 


0 


H represents column vector containing 
displacements and is given below 


0 


(3.8c) 


unknown nodal 


and 


W’- t 


u r U 2* U 3' *1' *2 ) 


}' 


(3.8d) 


The matrix [ B J® does not involve the nodal degree of 
freedom. Therefore the first variation of the linear membrane 


strain can be written as 



3C- 


*{ « } - c b if e “ * i 


3. 3. 3. Expression for d 


{-r 

- v Jnl 


Rewriting the Eqn. (3.7b) 


{«}“ 

i Jnl 


1 

2 


du 


3 0 




*u. 


dx. 


■1 

r > 


Su 3 


**■> 

H 

au 3 


dx 2 

J 

» - 


C3 


^3 ^3 

dx 2 Sx^ 

This nonlinear mid plane strain can be written in terms of 
dof as. 


r ^o 1 e e f y 

{-L ■ ;[ £] tG] W 


(3 


Where 


and 


[G] e = [G,, G, 


G 9 3 


[ G J 


m 


o 


Sx 


&H 


dx. 


(3 


(3 


(3.9) 


. lOa ) 


nodal 


. 10b) 


. lOc ) 


- lOd) 


The equation (3.10b) can be rewritten as 



H oi e 

■ - t b i:, t * i 

ML 2 


(3.11a) 


where [ B ] 
and 


o 

ML 


C B 01 r B 02 


r b ] 

i. OJ-* 


OJ-< ML 


0 


o 


* u 3 

f sn * 1 

J 

Sk ^ 
V 

t V 

3 ] 

.1 

. * x 2 . 

<?x 

^ cJ 


B Q9 3 ml 


0 0 


' * u 3 

dN 6 ! 

J 

4. 

*“3] 


V 




x . 
1J 


(3. 11b> 


0 0 


0 


0 

(3. 11c J 


The first variation of nonlinear mid plane strain can be 
written as 


1 [ [ b c { - •* } e + d c b ]^ e { * y ] <3 - i2 > 

e 

In [ B matrix, the term u 3 is expressed in terms of shape 

functions. After applying variation and rearranging the terms we 
obtain as, 

[ -[»£{*}■] - [ i *&{**}' ] (3 - 13 ' 


Substituting the relation <3*13) in to (3* IS) gives 
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r 

* t 

' e 


{ * } 
v. 1 ML 

C B ll { 

dd j 

(3.14) 


3.3.3 Expression for First, Variation of' Total Mid Plane Strain s 


The first variation of total mid plane strain can be written 
using relations (3.9) and (3.14) as. 


Where 


„{«}“ - [»]••{**}• 

[»]' - [ e»i:‘* [-]»• ] 


(3.15a) 


(3.15b) 


The total mid plane strain is given by. 


{-r ■ i {-y 


(3. 16) 


3. 3. 4 Expression for d | e J and d e | : 


The bending and shear strain at the mid plane of the 
laminate can be written again in terms of nodal displacements as. 


{-} = [c]e W 
{-}’ ■ t D]e W 


(3.17a) 
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Where 


and 
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O 




[ e] e 
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[ <v C 2 

— S ] 




(3.17b) 
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0 j 

0 1 


<9N 


“o 0x, 
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O S K 
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J 

Jx, 


«*N 


Sx 


0H 


&x. 


(3.17c) 


C D J 


N' 

R, 


N 

J 

R, 


0 


0N 

J 

x. 


«»IM 

J 

^x, 


N 


N 


(3- 17d) 


The first variation of bending strain and shear strain ii 
given by. 


i-r - 

[ C 

r 

£ d <5 

and *( . } s = 

[ D 

r 

d S 

Substituting variation 

of 

strain 


(3.18> 


(3.18) in to virtual work equation (3.5) we get. 
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r T 

N/] 

T j 

tO 

y IS ) = 


r -\ e t 

{dsj [ B] e - 

1 

«2 

h + { d<5 } E c l 6 ' 

1 

A J 


- 

N 6. 

V. J 

- »6 


+ 


{ 



dx^ dXg 


(3.19) 


Where R contains nodal forces and is 
3.19 can be written after taking out 


derived 



in next section. 


Eq. 



f N .1 



r r * 

1 

T 

1 

»(«■ [ B ] ■ 

N 2 

. + [C]«- 

V 

J 1 



M 6> 


0 

(3.20) 



Eq. (3.20) represents a set of non linear equations to be 
solved for a bending problem. 

3.3.5 Expressions For Stress Resultants s 


The expressions for stress resultants N^, Ng, N^» r Mg, 
M 6 , and Gg can be obtained by substituting the relations (3.16) 
and (3.17a) into Eq«(2.7) we get. 
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{ N 1’ N 2' N & } 

{ "l' M 2' M 6 } T - 

{v^r - 


[ ° fl ] 
t ° B ] 

[ aB J 
[ aD ] 

[ as] 


[c r{*y 

e<- -v e 1 
[ c ] e {' 5 } e 
c»]'w 


[ aA 3 [ B izl { * } e * 

(3.21a) 

[ 0B ] [ B “J { * } e+ 

(3.21b) 


(3.21c) 


Substituting the expressions 


for stress resultants ie.» 


relations (3.21a) , (3.21b) and (3.21c) in (3.20) we obtain 




J[ 


e T 1 


[ B ]° C SA] [ B ]“ + - [ ■>]*" [ «3fi] [ B]°* * 


[ B ]f [ •»] [C]' - [ Bit’ [ «A] [ B ]f + 


( [ ®] [ B C. * ( B 3l [ *] [ C]“ * 

t 1 e T e 

[ C ] e [ QB ] C B ]“' + - [ c ] e [ QB ] [ B ]“* + 


[ C ] [ aD ] [ C ] 


* [ d] bT [ as] [ D] e ] 

{ 5 } dx 1 d *z ■ { R } = ( 


(3.22) 
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3 . 4 . Solution of NonLinear Equilibrium Equations s 


The nonlinear equilibrium equations can be solved by two 
different methods viz. the direct iteration method with 
incremental approach and the incremental newton— raphson method. 

In the di rect iteration method load is applied in small 
increments. Some initial solution is guessed (typically, 
<5 = -|o|) and the equilibrium equations are solved to get new 
solution . This solution is used again to solve the equations . This 
type of iteration is continued till the two solutions converge . 
This method may not converge for highly nonlinear problems . 
Moreover it is computationally expensive. 

In the Newton Raphson method to start with a solution is 
guessed. An improved solution < <5^ + d <5^ ) is obtained from 


<5.+ d.5. ) 



[ « ] T { d * } 



(3.23a) 


where K is the tangent stiffness matrix given by, 

T 


[ “I, 


' *V(<5. ) 

au s) e 


(3.23b) 


Eqn. (3.23a) gives nonlinear incremental equilibrium 
equations which gives the linearized approximation to the relation 
between residual forces and the incremental displacements. 

The new improved solution can be found as 



37 


<5. — & + d <5 

vt-i v i 


( 3 » 24 > 


To improve numerical stability the 
increments . 


load 



is 


applied 


in small 


3. 5 Tangent, Stiffness Matrix a 


The tangent stiffness matrix for a element can be written 

as , 


[ K], 
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Where w is the total degree of freedom of the element. 
An element K in Eq. (3.20) can be written as 

Ti,j 
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Where 


*R 


<9«5j 


= 0 


It is shown in the Appendix — II that the tangent stiffness 
matrix can be written as 
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(3 

.27) 

Using the 

relations 

(3. 

15b) 

in 

Eqn. (3.27) 

the tangent 


stiffness matrix can be written as 


[ *] r - [ K ], + [ K ]«. + C K L <3 - 28a) 


Where 


K 


K 


ML 


K 


Linear stiffness matrix 
Initial displacement matrix 
Initial stress matrix 
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K I = 

Jo 1 


T 

[ G] e 
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A J 
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N 2 1 


( 3 . 28d ) 


3. 6 Elemental Force Vector * 


The work done by external loads for the bending problem may 
be independent or combination of these load cases. 

(i) Concentrated loads R £ acting in the di rect ion of the 



corresponding nodal degree of freedom. 

( ii ) Uniformly distributed load on the surface . 


4 0 


The total external work done by these forces on an element 
can be written as 


W, 


ER 


W, 


ER 


{ * }' { *<=} + I { 3 T R o dS 

t-r’ w 


(3.29) 


R £ is vector containing nodal forces/couples at node a in the 
direction of that nodal degree of freedom. The vector R q can be 
written as 


and R is 

OJ 


each S 

j 


{ r o} T - {'or 

given by Q 

{ s } = { ^ ’ 
is given by / 0 



(3.30) 


(3.31) 


3.7 Geometric Stiffness Matrix for the Plate Problem * 


For the linear buckling problem the geometric stiffness 
matrix can be derived by considering the work done by m plane 
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•forces N^and at the mid plane of the laminate. The datum 
configuration is taken as ini t ial ly compressed configuration. The 
mid plane strains due to bending can be written as 



The work done by these forces is given by 


(3.32) 





* U 3 1 

dx dx 

**2 J 

(3.33) 


Eqn. (3.33) can be written as 
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* U 3 1 
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1 


Sx A 
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N, 


N, 


N, 


-j 

. \ 
du 3 


&x A 
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dUg 

j 

‘ ^>X_ - 


dx. dx„ 
1 2 


(3.34) 


£U_ dUg 

Expressing the terms and — in terms of shape 

functions and nodal displacements Eqn . ( 3 . 34 ) can be rewri tten as 
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1 T T 

- [*] e [G] e 

A 

In this expression ,the stress resultants and are 

constant over the domain and will be equal to applied forces- The 
in plane forces N^, N^and can be expressed as 

N^l = X ; Ng = X n 2 ; = X (3.36) 

Substituting the relation (3.36) in (3.35) and rewriting the 
expression for external work we obtain 

W 6 = ^ X [ <5 ] [K] g [6] S (3.37) 

Where is the geometric stiffness matrix 

The first variation of the work done by external in plane 
forces can be written by considering the symmetry of as, 

du| -X [ d< ]•’ [ k]g I s )* (3 - 38 ’ 

Now the equilibrium equations for a linear buckling problem can be 
written as, 

[-*]• [«]; r [<•*]• [•‘ist'i *- 0 (3 - 39> 
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or [ C •< ] “ -X [K ]■][*]■. 0 (3.40) 

Eqn . (3-40) represents typical eigen value problem. On 
solving, the lowest magnitude of eigen value gives the buckling 
load. The eigen vector represents the mode shape. 

3.8 Geometric Stiffness Matrix for the Shell Problem : 


The geometric stiffness matrix or the shell problem can be 
derived by considering the work due to pre buckling stresses (N^, 
Ng, and N & ) during the buckling. 

The work done expression can be written as 



2N 6 



(3.41 ) 


Expressing u^ , Ug and u^ in terms of shape functions and the total 

expression will become similarly as Eqn. (3.37) .This K^. matrix is 

substituted in equilibrium equations and solved for the eigen 

value problem. The lowest eigen value is the buckling load. 

The integrals involved in the expressions of y(<5), \ T 

and K are evaluated using gauss - quardrature. The exact 
a 
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integration of the stiffness matrix requires 3x3 gauss point rule. 
But, at the thin plate limit to overcome locking 2x2 gauss 
quadrature is used for shear energy terms, while 3x3 is used for 
non shear terms . In the above section equilibrium equations are 
developed for an element. These equations can be assembled using 
the connectivity at the element interface. Proper boundary 
conditions are applied and the assembled equations are solved. 



CHAPTER 4 


RESULTS AND DISCUSSION 


4.1 Introduction s 


The equations developed in Chapter 3 are employed for 
numerical studies on some typical plates and shells. The results 
obtained ( ie. , present linear theory (PLT) , present nonlinear 
theory ( PNLT ) ) , are compared with some experimental , higher order 
theories ( HSDT ) and classical plate theory (CPT) results. Results 
are presented in the form of graphs and tables. The parameter 
whose effect is examined include number of layers, side to 
thickness ratio ( radius to thickness ratio for shell > and 
orthotropy. The effect of different simply-supported boundary 
conditions on bending-stretching coupling is also investigated. 

The finite element used for the discretization of plates 
and shells is a nine noded shell element. In the present work a 
shear correction factor ( 5 / 6 ) has been used. Reduced integration 
has been used for transverse terms in the case of thin plates and 
shells because full integration leads to over stiff results. The 
plate problems are solved by taking the radius of curvature very 
large (typically 10 30 ) in strain displacement relations Eqn. 
(2.3a) The boundary conditions used in the present work are shown 

in the Fig. 2.3. 
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For linear problems the matrices K and K t in formulating 
tangent stiffness matrix are put as zero. For stability problems 
the matrix K^_ matrix is formed and Eqn. (3.40) is solved for the 


eigen value. 

For nonlinear bending problems a solution <5 
assumed and for the first iteration linear solution 


= { ° } is 

is obtained. 


Using the 1 inear solution the matrices K and K are computed and 
solution is found again from the Newton - Raphson iteration. The 
iterations are continued till the two successive solutions 6 1 , 

j[ + •'I 

S converge. The following convergence criterion is used. 



Where i is the iteration counter and k is the nodal 
variable counter. 

Whenever symmetry exists quarter plate model has been 
employed. Unless specified, a 3X3 mesh is used for bending and a 
2X2 is used for buckling problems. The following nondimensional 
deflections and buckling loads are used. 
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Where u^ is calculated at the center of the plate> N is 
the applied critical load in case of linear bucfcling problem * 

4.2 Presentation of Results and Discussion s 


The following types of problems are solved using the 
formulation given in chapter 3. 

4.2. 1 Bending Analysis 

(a) Numerical accuracy of the present model with respect to 
the element 

( b) Bending analysis of three layered cross-ply 

( c) Bending analysis of angle-ply 

(d) A comment on the boundary conditions 

(e) Bending of isotropic spherical shel 1 

( f ) Bending of cylindrical shell 

( g ) Cross— ply freely supported cyl indri cal shell under 
sinusoidal load 

4.2.2 Stability Analysis 

(a) Stability of isotropic plate 

(b) Stability of symmetric cross-ply 

(c) Stability of angle-ply 

(d) Stability of isotropic cylindrical shell 

4.5.3 Nonlinear Bending Analysis of Plates 

Ca) Nonlinear analysis of orthotropic plate 

(b) Nonlinear analysis of symmetric cross-ply 

(c) Nonlinear analysis of simply supported cross-ply and 


angle— ply 
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(d) Nonlinear analysis of clamped corss-ply and angle-ply 

(e) Sandwich plates 

4.2.4 Post Buckling of Isotropic Plate 


4.2.1a Numerical Accuracy of the Present Model with respect to 


the Element : 


To check the numerical accuracy of the present element a 
test problem is solved for which closed form solution (CFS) is 
available C113. 

A simply supported square cross-ply (0/90/90/0) laminate 
with the following properties was considered under sinusoidal 
loading . 


= 25 


v 


12 


0.25 



= 0.5 



= 0.2 , 



(4.1) 


The results of non dimensional ized deflections u 3 were 
obtained for the following cases 

(a) Full integration (3X3 gauss rule for bending terms and 
3X3 gauss rule for shear terms) 

(b) Mixed integration (3X3 gauss rule for bending terms and 
2X2 gauss rule for shear terms) 

(c) Reduced integration (2X2 gauss rule for bending and 

shear terms) 

The solutions were obtained using si boundary conditions 
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which are used in Ref CUD for the quarter plate model. 

The results for various thickness ratios were obtained and 
presented in table 4.1. We can conclude from the results that 

(1) The nine noded element give virtually same results for full , 
mixed and reduced integration schemes . 

(2) As thickness of the plate is reducing it can be seen that full 
integration gives over stiff results. This is due to the reason 
that the present model being not able to satisfy the Kirchoff *s 
condition as a constraint. It can also be seen that the results 
with mixed integration scheme exactly agree with closed form 
solutions . No signif i cant dif ference is observed between reduced 
integration and closed form solutions. However the experience of 
the people in this field says that mixed integration is suitable 
for thin plate or thin shell situations. Therefore we used reduced 
integration for the shear terms in these cases. 

(3) The finite element model used here is expected to yield 
correct results. 

4.2.1b Bending analysis of Three Layered Cross-Ply s 


The object of this study is to compare the present theory 
with the closed form solutions of higher order shear deformation 
theory (HSDT) presented by Reddy C163. The higher order theories 
are computationally expensive and involve large number of degrees 
of freedom. 

A simply supported three layered corss-ply (0/90/01 of the 
following material propert ies subjected to uniformly distributed 



lateral load is taken up for which HSDT results are available 

E 1 = 19.2 x 10 6 psi T Eg = 1.56 x IQ 6 psi, 

^12 = ^13 = x 1^ P s i i Ggg = 0.523 x lO^psi, 

v ^2. ~ » a = b = 10 in. (4.2) 

The results are obtained using si boundry conditions. The 
non dimensional i zed def lections u^ are presented for various 
thickness ratios in table 4.2. 

It can be seen from the results presented in table 4.2 that 
the present results are close to the HSDT results for a/h >10 
without much computational effort . The maximum difference is about 
2.7%. The present results underestimate the deflections compared 
to HSDT. It can be seen that the difference reduces as the 
thickness of the plate decreases. The CPT results underestimate 
the deflections and are not reliable at thick plate limit. It can 
be concluded that present theory is adequate in predicting the 
deflections for a/h > 10 

4.2.1c Bending Analysis at Angle-Ply : 


In this section the effect of bending— stretching coupl ing 
on linear deflections for an angle-ply (-45/45/...) is studied. A 
simply supported angle— ply made of graphite epoxy of the following 
material and geometrical properties subjected to uniformly 
distributed lateral load is taken up. 
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= 40 , G 23 = 0.5 E 2 , 

= 0.25, E g = 10 6 psi 

b =10 in. 


12 


G 13 0,6 E 2* 


(4.3) 


Fig . 4.1 shows the plot of nondi mens i ona 1 i zed central 
deflections for various thickness ratios with increasing number of 
layers . The results are obtained using s2 boundary conditions in a 
quarter plate . It can be seen that the deflections for two layers 
are higher than the four layer and eight layers . This is due to 
the bending stretching coupl ing which is maximum for two layers 
and decreases as the number of layers is increased. 

4. 2. Id A Comment, on the Boundary Conditions : 


To study the effect of different in plane boundary 
conditions on bending— stretching problems a two layer angle— ply 
and cross— ply of the material properties given in Eqn. (4.3) 
subjected to uniformly distributed lateral load is taken. 

a = b = 10 in. h = 0.1 in. C4.4) 

One can use si boundary conditions and s2 boundary 
conditions as a symmetry conditions for the quarter plate model 
associated with cross— ply and angle— ply respectively. However one 
cannot use the symmetry conditions si in the quarter model of 
angle-ply and s2 in the quarter plate model of cross-ply. The 
effect of different boundary conditions on the linear bending of 

CENTRAL L'RRARY 

' ' t.. ; xA-vruf* ■ : 



40 , 


G 23 = °' 5 E z 



r 


0.25, Eg = 10 6 psi 


a = fa = 10 in. 


G „ = G,„ 0.6 E_, 
12 13 2 


(4.3) 


Fig. 4.1 shows the plot of nondimensionalized central 
deflections for various thickness ratios with increasing number of 
layers . The results are obtained using s2 boundary conditions in a 
quarter plate . It can be seen that the deflections for two layers 
are higher than the four layer and eight layers. This is due to 
the bending stretching coupling which is maximum for two layers 
and decreases as the number of layers is increased. 


4. 2 . Id A Comment on the Boundary Conditions s 


To study the effect of different in plane boundary 
conditions on bending-stretching problems a two layer angle-ply 
and cross-ply of the material properties given in Eqn. (4.3) 
subjected to uniformly distributed lateral load is taken. 


a = b = 10 in. h = 0.1 in. 


(4.4) 


One can use si boundary conditions and s2 boundary 
conditions as a symmetry conditions for the quarter plate model 
associated with cross-ply and angle-ply respectively. However one 
cannot use the symmetry conditions si in the quarter model of 
angle-ply and s2 in the quarter plate model of cross-ply. The 
effect of different boundary conditions on the linear bending of 
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cross-ply and angls-ply are showed in the table 4.3. It can be 

seen that for N = 2 layers, th* rfiffar.-, , * 

layers , ihe difference between the deflections 

obtained from these boundary conditions (SI and S2) is maximum. 

But for N = 10 the difference is small. This due to the fact that 

as the number of layers are increased the bending-stretching 

coupling decreases and the effect of different in plane boundary 

conditions is not significant. 

4.. 2. ie Bending of Isotropic Spherical Shell s 


A simply supported isotropic spherical shell shown in 
Fig. (4.2a) of the following material and geometrical properties 
is studied under a concentrated load of 100 lbs. 

7 

E = 10 psi, v = 0.3, h = 0.1 in. R 1 = R g - 96 in. 

a ss b = 32 in. (4.5) 

Due to the existence of the symmetry a 2X2 mesh is used to 
discretise one quadrant of the shell. As the shell is thin (R/h = 
320), mixed integration is used. The linear deflection obtained 
under the concentrated load is 0.03723 in. The mesh is refined to 
3X3 and the central deflection obtained is 0.03886 in. in 
comparison to 0.038661 in. of Rao C371 whose result is based on 
classical shell theory where shear deformation is neglected. 

4. 2. If Bending of Cylindrical Shell : 

A clamped isotrpic cylindrical shell shown in Fig. (4.2b) 
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of the following material and geometric properties subjected to 
uniformly distributed lateral load is considered for the study. 

E = 3103 N/mm 2 , v = 0.3, h = 3.175, a = 254 mm. 

R 0 = 1000 * IQ" 6 N/mm 2 & = 0.1 rad. (4.6) 

Here a 3X3 mesh is used to descretize one quadrant of the 
shell. Clamped boundary conditions given in Eqn. (2.10) are used. 
The linear central deflection obtained at the center of the shell 
is 0.S4297 in. compared to 0.8250 in. of Dutt whose result is 
based on classical shell theory. The Dutt result was given in 
Ref .C393 


4. 2. ig Cross-ply Freely supported Cylindrical Shell under 

Sinusoidal Load i 

A cylindrical cross-ply (0/90/90/0) cylindrical shell of 
the material properties given in Eqn (4.1) subjected to 
sinusoidal loading is taken for the study. The geometrical 
properties are 

30 

a/b = 1.0, a/h = 10 R g “ W 

fix* nx„ 

R * B SIN — 1 SIN -r (4 - 7) 

o o a a 


The following boundary conditions are imposed 



54 


U 3 = u 2 = 0 ** ^ = 0,a 

u 3 = U 1 * 0 at = 0,b ( 4 . 8 ) 

For this special case of boundary conditions , loading closed form 
solutions are available in Ref £382 with the same theory. A 4x4 
element mesh is used to descretize the quarter shell. The 
nondimens ional deflections u^ are presented in Table 4.4 for a/h = 
10 and a/h —100. It can be seen that the FEM results are close to 
the closed form solutions. 

4. 2. 2 Stability Analysis : 

In this section stabi 1 i ty of isotropic plates and cross— ply 
plates have been compared with HSDT results. Bending stretching 
coupl ing on angle-ply laminates has been studied. Results are 
presented for cylindrical shell under axial compressive load. 

4,. 2. 2a Stability of Isotropic Plates : 

The object of this study is to evaluate higher order 
theories with respect to their computational power and accuracy- A 
simply supported isotropic square plate of the following material 
and geometric properties under uniaxial compressive load is 
considered. 

E *= 10 7 psi, v = 0.3, a = b = 10 in. (4.9> 

A 2X2 mesh is used to discretize the quadrant of the plate, 
si boundary conditions are used and linear buckling load is 



computed by solving the ei gen v - ri ._ „ 

alue problem. Nondimensional 

buckling loads are presented i n table 4 5 , 

form solutions of HSDT by Reddy C17H for various thickness ratios. 
It can be seen from the results the present results are very close 
to HSDT results for all thickness ratios. The maximum difference 
is 0.75%. The difference is decreasing as the thickness is 
decreasing. But in the results shown this effect cannot be seen 
over the whole range due to the small error in the eigen value 
calculation of large matrices and locking phenomenon. The CRT 
results over estimate the buckling loads for thick plates (a/h < 
20) . This due to the shear deformation neglected in CRT. 


4 . 2. 2b Stability of Symmetric Cross -Ply i 


A simply supported graphite epoxy symmetric cross-ply 
(0/90/90/0) laminate of the material properties given in Eqn. 
(4.3) is considered to study the uniaxial buckling load. 

The nondimensional uniaxial buckling load N x was computed 
and presented in table 4.6 with HSDT results C173. si boundary 
conditions are used. 

It can also seen from the results that the present results 
are very close to HSDT results. The maximum difference is 1.5%. It 
implies the adequacy of the present theory in predicting the 
buckling loads. HSDT in FEM will involve a lot of additional 
computational effort without improvement m the results. From 
Table 4.5 and 4.6 we can conclude that the difference between the 
present theory and the HSDT is slightly increasing 



orthotropy. 


4. 2. 2. c Stability of Angle-Ply : 


A simply supported square angle-ply (45/-45/. » . ) under 
uniaxial compression load is considered for the study. The 
properties of individual laminate are given in Eqn. (4.3> m 

Solutions are obtained for S2 boundary conditions in a 
quarter plate. Fig. 4.3 displays the variation of nondimens i anal 
buckling load 13^ with fibre orientation and number of layers. It 
is shown C413 that the solution approaches to orthotropi c solution 
at N = 6. Hence in the present study N = 6 and N = 2 are 
considered for a/h = 10 and a/h = 100. It can be seen that as the 
number plates increased the buckling load is increasing for the 
same thickness. This is due to the reduction in the bending 
stretching coupling and buckling load is minimum for two layer 
laminates. There is significant difference between the results for 
a/h =* 10 and a/h =* 100. This is due to shear deformation which 
decreases the buckling load. At 45° orientation the buckling loads 
are maximum for all number of layers. 

4. 2. 2. d Buckling of Isotropic Cylindrical Shell under Axial Load s 

An isotropic homogeneous cylindrical shell with length 11 
in., radius of 36 in., and thickness of 0.125 in. is considered as 
a test problem to check the formulation. Both the edges are simply 
supported. The analytical solution for the same is available with 
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Donnell * s theory in Ref C423. For the finite element modeling it 
is assumed that symmetric conditions prevails at the mid length 
and that the cylinder buckles in 15 circumferencial full waves as 
shown in the Fig. (4.4a). Consequently a 5.5 in. long and 12° wide 
cylindrical segment is considered with symmetrical conditions 
imposed on the three sides. The following boundary conditions are 
used 



= Ug 

= *2 

= 0 

DO 

x i 

li 

Q 

w 

fil 


u i 

3 

1! 


= 0 

on 

X 2 

= 0 

o 

V 

a 

w 

u z 

= * 2 

= 0 


on 

X 2 

= b 



where coordinate axes pass through left corner of the segment as 
shown in the Fig (4.4b) . 

The eigen value problem of Eqn. (3.42) is solved and 
buckling load is applied as 2721 lb/in. The result is with a 
difference of 3.5% with the analytical solution (2625 lb/in.) 
given in Ref .C423. This difference is likely to get reduced as the 
mesh is refined. This is being the reason that as the shell is 
very thin and the shear deformation effect is expected to be 
negligible and the results obtained from the present theory can be 
compared with the Donnell's theory given in Ref. E423 which does 
not take into shear deformation. 

In the present study we could not make use of a routine 
SSPACE given in Ref. C433 which computes specified number of 
lowest eigen values only. Hence we could not go in for a refined 
mesh than 2X2 because of the computational limitation. Laminated 
shells buckles in more than one half wave along the axial 



direction- Hence the buckling analysis laminated shells is not 
taken up. 


4.2.3 Nonlinear Bending of Composite Plates s 


IN this section the finite element results obtained from 
the present theory are compared with the experimental results and 
classical nonlinear results. As the bending stretching coupling is 
maximum for two layer laminates the influence of different 
boundary conditions is expected to be more pronounced. To 
investigate this the solution is obtained with SI, S2 and S3 
boundary conditions for two layer angle— ply and cross— ply. 


4.2. 3. a Nonlinear Analysis of Grthotropic Plate s 


A simply supported orthotropic plate of the following 
material properties subjected to uniformly distributed load is 
chosen for which experimental results and classical nonlinear 
results (CNLT) are available in Ref. C183 


«& = fa ^ 12 in « 

E. = 3 x iO 6 psi. 

1 

G 12 = G 23 = G 13 

1 v 33 U*' 255 ■ ’ 23 m 

v 12 23 13 


h = .138 in. 
E 2 = 1.28 x 10 6 
.37 x 10 6 

(4. ID 


4.2.3a CD 
The 


Newton- Raphso n vs Picard Iteration s 
solution of this problem with SI boundary 


conditions 



'-j o 

was obtained by using both Newton- Raphson and picard iteration. 
Tabie 4.7 shows the results obtained from Newton-Raphson and 
picard type of iteration. 

It is evident from Table 4.7 that the number of iterations 
for convergence increases with the load when picard iteration is 
used. However Newton-Raphson converges in three iterations for all 
loads. Hence in the problems investigated below Newton — Raphson 
method has been employed. 

A 2X2 mesh is used in quarter plate. Using Newton-Raphson 

results were obtained with SI and S3 boundary conditions. Fig. 

(4.5) displays the results obtained from present linear results 

(PLT) , present nonlinear results ( PNLT ) , classical linear results 

(CLT) and classical nonlinear results (CNLT) given in Ref . Cl 83. 

It can be seen that the PLT and CLT results are close. This is due 

to the fact that as the plate is very thin (a/h = 87) the effect 

of shear deformation is expected to be negligible. It can be seen 

that compared to CNLT results the PNLT results obtained with SI 

boundary conditions are in close agreement with the experimental 

results. This is due to the reason that the boundary conditions 

used in Ref. CIST to obtain the CNLT solution are u i “ u 2 

u_ = 0 and ML *= O at the simply supported edge which may not 
3 E 

correspond to the the experimental conditions. The PNLT results 
obtained with S3 boundary conditions are away from the 
experimental results. The difference between SI and S3 boundary 
conditions are due to the bending stretching introduced through 
the nonlinear term. 

It can be concluded that the practical simply-supported 



boundary conditions for the orthotropic plate can be in between SI 
and S3 boundary conditions. 

4 , 2. 3b Nonlinear analysis of Cross-Ply s 


A clamped symmetric cross-ply (0/90/90/0) of the following 
material and geometric properties used to study the nonlinear 
analysis under uniformly distributed load. For this data CLT, CNLT 
and experimental results are avalible in C183. 


b = 12 in, 


E„ = 1,828 x 10 
1 


h = 0.096 in. 

Eg = 1.8315 x W 




G ,2 = G 13 " G 23 = 3 ' 1£5 * 10 
« " V 13 =“23 - ■ 2395 


< 4 . 12 ) 


The CPT, CNLT, PLT, PNLT end experimental results are 
presented in Fig. (4.4» with load intensity. In this case there is 

a difference between CLT and PLT reults even though the plate is 

. - th _ n-f C183 the following boundary 

thin. This is because m the Ref. 

. * th _ riamoed edge in the linear and 

conditions are used at the clamp y 

nonlinear analysis 


du. 


(4.13) 


Ui = u 2 = u 3 


ax. 


which over estimates the deflections. But in the present FEM model 
it is possible to specify the correct boundary conditions which 

• p„„ (2 10) The PNLT results are lower than the 

are given in Eqn ■ U- • 

. , trend. The reason of 

experimental results which shows e 
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difference in the results is that the clamped boundary conditions 
are not satisfactorily imposed during the experimentation as 
stated in the Ref .£183. The CNLT results shows altogether 
different trend because of the reason mentioned above regarding 
boundary condi t i ons . These two examples led to the confidence in 
the present model in predicting the accurate nonlinear response of 
the composite plates . 

2- 3c Nonlinear Analysis of Simply Supported Cross— Ply and 
Angle-Ply s 


A simply supported square plate of the following material 
properties is taken for the study. 



„ Eg = 7.031 x 10 5 N/cm 2 , 


23 


= .2 


V 12 


243.8 cm. 


h = .635 



.5 


.25 


(4.15) 


As the bending— stretching coupling is maximum for two layer 
laminates a two layer angle-ply and cross-ply subjected to 
uniformly distributed lateral load is considered. The load 
deflections curves are shown in Fig. (4.7) and Fig. (4.8). For the 
cross-ply SI and S3 boundary conditions are used. It can be seen 
that S3 boundary conditions are stiffening type than SI boundary 
conditions. The same trend is shown for the orthotropic plate in 
Fig. (4.4). For the angle-ply S2 and S3 boundary conditions are 
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used- Here also S3 boundary conditions show the same trend with S2 
boundary condition. 


-4. 4,. 3 d Nonlinear Analysis of Clamped Cross-ply and Angle-Ply s 


A c lamped cross— ply and angle— ply of the material 
properties shown in Eqn. (4.3) is used to study the nonlinear 
behavior . 

Fig . (4.9) and Fig. (4. 10) shows the load deflection curves 
for two layer cross-ply and angle-ply. It can be seen that the two 
layer cross-ply exhibits more nonlinearity than angle-ply. It can 
also seen that for cross-ply as the number of layers are increased 
from two layers to eight layers the nonlinearity is decreasing. 
This is due to the reduction in bending-stretching coupling with 
the number of layers . 

4. 2. 3© Sandwich Plates s 


The object of this study is to validate the present theory 
to sandwich plates by comparing with the available experimental 
results. A simply supported square sandwich laminate with 
aluminium faces and balsa core subjected to uniform transverse 
pressure is studied with the following face properties. 

E f = 10 x 10 6 , t f = 0.32 in. 

A 2x2 mesh is used to idealize one quarter plate and S3 
boundary condition are used. The results are obtained from the 
present theory and experimental results are presented in Table 4.8 



and table 4.9 for varying spans, core thickness, core shear 
rigidity and load intensities. The present nonlinear results 
compare satisfactorily with the experimental results. When the 
shear rigidity of the core high it can be concluded that 
neglecting the face shear rigidity yields better results than 
considering the shear rigidity of the face plate. 


4. £« 4. Post buckling of Isotropic Plate : 


A simply supported isotropic square plate of the following 
material properties is chosen for the study 


E = 30 x 10 6 Psi, a = b = 10 in. h = 0.1 in. 

In this problem following type of boundary conditions were 


used 




nx 


« Kg 


= c 1 

SIN 

I 

cos 





nx 4 


rcxg 


= C 2 

COS 

1 

2a 

SIN 

b * X 2 

(4.16) 





nx. 


" C 3 

COS 


COS 

-ar 



Here the origin is chosen at the centre of the plate. 2a and 2b 
are the dimensions of the plate respectively. 

The boundary conditions can be rewritten as 


Ug = Xg. 

u 2 = + « b 


on all edges, 
on X/j = i a 
on x„ = 


(4.17) 


+ b 



After the application 


of these boundary conditions the 
problem becomes the minimization of strain energy as the variation 
of potential energy is zero. 

From Ref. C423 the critical value of € where a buckled 
configuration exists for v = 0.3 is 


■fi = 0 . 632 — = 

a^ 


<4. IS) 


ie., if « < < 2 ^ ^the plate will be in compressed configuration only 


Let n = 

<6 

cr 

A small lateral load 21bs (0.18% of the critical load) is 
appl ied at the centre of the plate to initiate buckling. The 
buckled configuration serves as starting vector for the iteration. 

The central deflection computed for various values of n 
iteratively are presented with analytical Gelerkian solution 
obtained in Ref. C423 in Fig. 4.11. It can be seen that the 
present FEM results and Gelerkian solution are agreeing. 















Non dimen s ipnal ized deflections for symmetric cross-ply (0/90/90/0) 
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r^-Lle 4-‘7 Results obtained for Newton-Raphson and 
and Picard Iteration for the Orthotropic 
Plate 


Load 

No . of 

iterations 

Picard 

Newton-Raphson 


1 

4.5365 E-02 

4.5365 E-02 


2 

4.3818 E-02 

4.4142 E-02 

0.2 





3 

4.3928 E-02 

4.3917 E-02 


4 


4.3928 E-02 


1 

8.7836 E-02 

8.5241 E-02 


2 

8.0171 E-02 

8.2001 E-02 

0.4 

3 

8.1796 E-02 

8.1515 E-02 


4 

8.1467 E-02 

8.1522 E-02 


5 

8.1578 E-02 

8.1522 E-02 


1 

0.1223 

0.1154 


2 

0. 1085 

0.1125 

0.6 

3 

0.1135 

0.1122 


4 

0.1117 



5 

0.1123 



6 

0.1121 



1 

0.1497 

0.1399 


2 

0.1316 

0.1377 


3 

0.1405 

0.1376 


4 

0.1362 


0.8 

5 

0.1383 



6 

0.1372 



7 

0.1377 



8 

0.1375 



1 

0.1720 

0.1608 


2 

0.1514 

0.1592 


3 

0.1638 

0.1591 


4 

0.1563 



5 

0.1607 


1.0 

6 

0.1581 



7 

0.1597 



8 

0.1593 




0.1589 



9 











Table 4.8 Sandwich plates comparison with experimental results ( G c =12670ps i ) 
A - G.not considered B - G.cons i dered 
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Table 4.9 Sandwich plates comparison with experimental results (G =5840psi) 
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CHAPTER 5 


CONCLUSIONS 


5.1 Conclusions : 

A first order shear deformation theory for the nonlinear 
bending , post buckling of plates, linear bending and buckling of 
shells is presented. On the basis of the limited study carried 
pout the following conclusions can be drawn. 

1) It is clear from the Table 4.1 that the type of 
integration scheme has a little effect on the nine noded 
lagrangian element . It is also evident that one has to use reduced 
integration for the shear terms in the thin plate region to 
circumvent locking . 

2) It can be seen from the Table 4.E that FSDT results are 
not rel iable for a/h < 10 for the bending problems. For a/h > 10 
there is no significant difference between FSDT and HSDT results . 

3) The FSDT results for the uniaxial buckling load are 
studied for three layered symmetric cross-ply and is given in 
Table 4.6. For this type of laminate studied it can be seen that 
present results are very c lose to HSDT results for all thickness 
ratios. It can be concluded that the present theory is adequate to 
determine the lowest buckl ing load. The classical plate theory is 
not reliable for a/h < 20. 

4> it can be seen from the Table 4.7 that Newton-Raphson 
method is more efficient than the pi card type of iteration through 
the former involves computation of derivative of the stiffness 
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matrix. 

5) As seen from the Fig. 4.4 it can be concluded that the 
present nonlinear results are close to the experimental results 
with SI boundary condition than S3 boundary condition . The 
difference in results for SI and S3 boundary conditions is due to 
the bending-stretching coupl ing introduced through the 
nonlinearity. The linear results over estimate the deflections and 
are not rel iable where the deflections are of the order of the 
thickness. It can noted that the experimental boundary conditions 
may lie in between SI and S3 boundary conditions for an 
orthotropic plate. 

6) Due the bending stretching coupl ing the different in 
plane boundary conditions (SI ,S2,S3) give different results . Of 
al 1 the boundary conditions S3 is stiffening type. This can be 
seen from Fig. 4.7 and 4.S. 

7) From Fig . 4.9 and 4.10 it can seen that cross— ply 
exhibits higher nonlinearity than angle-ply laminates. It can also 
noted that for the cross— ply as the number of layers is increased 
the nonlinearity is decreasing. 

5. 2 Recammendati ons for the Future Work : 

Some of the future work that can be carried out as an 
extension or improvement over the present investigation are 
suggested below. 

The buckling formulation and program that is developed can 
be utilized to analyze the stability of sandwich plates. The 
program developed here also caters for biaxial and shear buckling 
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of laminated composite plates. This could he investigated for 
future work. 

With the little effort this work can be extended to 
nonlinear analysis by taking al 1 the nonlinear terms in the strain 
displacement relation as 
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There will be additional matrices such as or give in 
Eqn « (3. ) . This becomes computat ional ly expensive. One can 
investigate whether that effort is required in engineering design. 

The post buckling problem done in the present work is 
1 imited to some special boundary conditions. By adding an 
additional external potential term one can solve it for al 1 
boundary conditions . 
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APPENDIX I 


1 Introduction s 


A composite laminate is composed by bonding together a 
number of lamina. In each lamina material axes may not coincide 
with the reference axes . The material constants of a lamina which 
are functions of engineering constants and orientation of fibre « 
are given in the following section. 

2 Description s 

The stress strain relations for a lamina, with material 
axes coinciding with reference axes and assuming the transverse 
stress is negligible, are given by 
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For an arbitrarily oriented laminated lamina after carrying 
out necessary transformation stress-strain relations are given hy T 
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In which C= COS (a) and S = SIN(oO and « is the angle between 
reference axes and principal orthot ropi c axes 



APPENDIX II 


1. Derivation of Tangent Stiffness Matrix : 


Substituting the Eqn . (3.20) , in to Eqn. (3.25) enables 
the tangent stiffness matrix to be written as 
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Where n denotes total degree of freedom. 

Since the stress resultants are 1 inear function of 
the variation of stress resultants can be writen 
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Eqn.<2.7) as 
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Substituting the relation (3. 15a) * (3. 15b) we can rewrite 
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Substituting the Eqn. (A 2.7) in to CA 2.1) we obtain the 
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expression for the tangent stiffness matrix as 
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Considering the third integral in Eqn. (A 2.8) and noting 
that [ B ]° e is not a function of <5 gives, upon mult i pi cat ion of 
l“l£ N ) and E 
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The differential with respect to S in Eqn . (A 2.13) are the 
coefficients of [ G ] matrix of Eqn. (3. lOd) typically for a 
degree of freedom l 
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Consequently Eqn . (A 2.17) may be reaaranged to give the 


expressions for K shown in Eqn C3.27a) 
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